"""
LSD 直线段检测算法
输入：灰度图像
输出：

"""
import cv2
import numpy as np
import math
import copy
import misc_func as mf
import LSD


class Rect:
    def __init__(self):
        self.x1=0; self.y1=0    #矩形主轴反方向上的端点
        self.x2=0; self.y2=0    #矩形主轴正方向上的端点
        self.width=0        #宽度
        self.cpx=(0,0)      #矩形的中心点
        self.theta=0        #矩形的方向（弧度）
        self.precision=0
        self.p=0
        self.dx=0; self.dy=0

class LSD_x:
    def __init__(self):
        self.NOTDEF=-1024
        self.M_3_2_PI = 3*np.pi/2
        self.M_2__PI = 2*np.pi
        self.scale=0.8  #图像的缩放比例
        self.sigma_scale=0.6
        self.n_bins=1024    #梯度伪序的bin数量
        self.ang_thresh= 22.5   #角度阀值
        self.log_eps=0.0          #阀值: -log(NFA)>log_eps
        self.quant=2.0             #受限于梯度范数上的量化误差
        self.density_thresh=0.7     #判断矩形对齐点密度的阀值
        return


    def log_gamma(self,x):
        """
        计算对数的gamma函数
        """
        q=[ 75122.6331530, 80916.6278952, 36308.2951477,
            8687.24529705, 1168.92649479, 83.8676043424,
            2.50662827511]
        if x>15.0:
            res=0.918938533204673 + (x-0.5)*np.log(x)-x \
                +0.5*x*np.log(x*np.sinh(1/x)+1/(810.0*x**6) )
        else:
            a=(x+0.5)*np.log(x+5.5)-(x+5.5);
            b=0
            for n in range(7):
                a=a-np.log(x+n)
                b=b+ q[n]*x**n
            res=a+np.log(b)
        return res

    def angle_diff(self,a,b):
        a=a-b
        while a<=-np.pi: a=a+self.M_2__PI
        while a>np.pi: a=a-self.M_2__PI
        return a

    def sub_sample(self,img, scale=0.8):
        """
         对图像进行减采样
        """
        if scale>1.0: scale=1.0
        height,width=img.shape[:2]
        new_width=np.int16(np.round(width*scale))
        new_height=np.int16(np.round(height*scale))
        out=cv2.resize(img,(new_width,new_height),interpolation=cv2.INTER_LINEAR)
        return out

    def gaussian_kernel(self,size,sigma,mean):
        sum = 0.0
        kernel = np.zeros(size)
        for i in range(size):
            val = (i - mean) / sigma
            kernel[i] = np.exp(-0.5 * val ** 2)
            sum = sum + kernel[i]
        if (sum >= 0.0):
            for i in range(size):
                kernel[i] = kernel[i] / sum
        return kernel


    def gaussian_sampler(self,image,scale,sigma_scale):
        ysize,xsize=image.shape[:2]
        N=np.int(np.ceil(xsize*scale))
        M=np.int(np.ceil(ysize*scale))
        aux=np.zeros([ysize,N])
        out=np.zeros([M,N])
        if scale<1.0:
            sigma=sigma_scale/scale
        else:
            sigma=sigma_scale
        prec=3.0
        h=np.int(np.ceil(sigma*np.sqrt(2.0*prec*np.log(10)))  )
        n=1+2*np.int(h);    #高斯核的尺寸

        #####先对X轴减采样
        double_x_size=2*xsize
        double_y_size=2*ysize
        for x in range(N):      #x为新图像中的x坐标
            xx=x/scale  #与x对应的原图像中的坐标
            xc=np.int(np.floor(xx+0.5))     #对xx的四舍五入,核操作的中心像素
            ##生成高斯核
            mean= h+xx-xc;
            kernel=self.gaussian_kernel(n,sigma,mean)
            for y in range(ysize):
                sum=0.0
                for i in range(n):  #以xc为中心的像素与核进行运算
                    j=xc-h+i;
                    while j<0: j=j+double_x_size
                    while j>=double_x_size: j=double_x_size-j
                    if j>=xsize: j=double_x_size-1-j
                    sum=sum+image[y,j]*kernel[i]
                aux[y,x]=sum
        ###经过aux之后得到ysize*N的图像，X轴方向减采样了
        ######对y轴减采样
        for y in range(M):
            yy=y/scale
            yc=np.int(np.floor(yy+0.5))
            ##生成高斯核
            mean = h + yy - yc;
            kernel = self.gaussian_kernel(n, sigma, mean)
            for x in range(N):
                sum=0.0
                for i in range(n):
                    j=yc-h+i
                    while j<0: j=j+double_y_size
                    while j>=double_y_size: j=double_y_size-j
                    if j>=ysize: j=double_y_size-1-j
                    sum=sum+aux[j,x]*kernel[i]
                out[y,x]=sum
        return out

    def levelline_angle(self,scaled_image,threshold, n_bins):
        """
         计算图像的梯度，等高线方向，以及进行梯度伪序操作
         threshold及论文中的rho
         返回值：梯度幅值，梯度方向，梯度伪序排序的像素坐标列表
        """
        ##计算梯度幅值，使用2x2窗口
        height,width=scaled_image.shape[:2]
        #out_image=np.zeros((height,width))  #输出图像
        mag=np.zeros((height,width))    #梯度幅值
        ll_angles=np.ones((height,width))*self.NOTDEF   #等高线方向的角度（弧度）
        max_grad=0.0
        ##计算梯度，角度
        for r in range(height-1):
            for c in range(width-1):
                com1=scaled_image[r+1,c+1]-scaled_image[r,c]
                com2=scaled_image[r,c+1]-scaled_image[r+1,c]
                gr=com1-com2;   gc=com1+com2
                #计算梯度幅值
                norm2=gr*gr+gc*gc
                norm=np.sqrt(norm2/4)
                mag[r,c]=norm
                #计算梯度角度
                if norm>threshold:
                    ll_angles[r,c]=np.arctan2(gc,-gr)
                    if norm>max_grad: max_grad=norm
                else:
                    ll_angles[r,c]=self.NOTDEF    #小于阀值的像素角度设为未定义

        ##得到梯度伪序
        #先创建列表
        bin_list=[]
        for i in range(n_bins):
            bin_list.append([])
        #按梯度值将像素分到对应的bin中
        for c in range(width-1):
            for r in range(height-1):
                norm=mag[r,c]
                i=np.uint(norm*n_bins/max_grad)     #计算像素梯度值对应的bin
                if i>=n_bins: i=n_bins-1
                bin_list[i].append((r,c))
        #按bin从高到低，将所有元素排成一列
        list_p=np.zeros([(width-1)*(height-1),2],dtype=np.int)   #
        n=0
        for i in range(n_bins-1,-1,-1):
            for px in bin_list[i]:
                list_p[n]=px
                n=n+1
        return mag,ll_angles,list_p

    def isaligned(self,x,y,angle,ll_angles,  prec):
        """
        判断(y,x)像素与给定区域角度值是否一致
        """
        a=ll_angles[y,x]
        if a== self.NOTDEF: return False
        angle=angle-a
        if angle<0: angle = -angle
        if angle>self.M_3_2_PI:
            angle=angle-self.M_2__PI
            if angle<0: angle=-angle
        return angle<=prec

    def region_grow(self,x,y,ll_angles,reg,used,prec):
        """
        区域生长算法，以(x,y)点开始进行区域生长
        ll_angles:等高线角度图，used:用于标记像素是否使用
        返回值：区域像素列表，区域角度
        同时会修改used
        """
        height,width=ll_angles.shape[:2]
        reg[0,1]=x; reg[0,0]=y
        px_count=1
        used[y,x]=1
        reg_angle=ll_angles[y,x]  #标量，记录整个区域的角度
        sum_sin=np.sin(reg_angle)
        sum_cos=np.cos(reg_angle)

        i=0
        while  i<px_count:
            #取一个元素进行区域操作
            cy,cx=reg[i]
            xx=cx-1
            while xx<=cx+1:
                yy=cy-1
                while yy<=cy+1:
                    if yy>=0 and yy<height and xx>=0 and xx<width and (not used[yy,xx]) and   \
                        self.isaligned(yy,xx,reg_angle,ll_angles,prec):
                        ##添加像素点到区域
                        used[yy,xx]=1;
                        reg[px_count,0]=yy
                        reg[px_count,1]=xx
                        px_count=px_count+1
                        ##更新区域的角度
                        sum_sin=sum_sin+np.sin(ll_angles[yy,xx])
                        sum_cos=sum_cos+np.cos(ll_angles[yy,xx])
                        reg_angle=np.arctan2(sum_sin,sum_cos)
                    yy=yy+1
                xx=xx+1
            i=i+1
        return reg_angle,px_count

    def region2rect(self,rect,region,px_count,reg_angle, precison, p, grad_mag):
        """
        为区域生成矩形
        reggion:区域像素列表  reg_angle:区域角度
        prec: 精度(弧度)    p: 像素是对齐点的概率，
        返回值：矩形的中心点，方向，宽度，长度
        """
        ##计算区域质心，作为矩形中心
        sum_x=0.0; sum_y=0.0; sum=0.0
        for i in  range(px_count):
            x=region[i,1]; y=region[i,0]
            weight=grad_mag[y,x]
            sum_x=sum_x+x*weight
            sum_y=sum_y+y*weight
            sum=sum+weight
        cx=sum_x/sum
        cy=sum_y/sum
        ##计算矩形方向
        Ixx=0.0; nIxy=0.0; Iyy=0.0; w_sum=0.0
        for i in  range(px_count):
            x=region[i,1]; y=region[i,0]
            weight=grad_mag[y,x]
            w_sum=w_sum+weight
            Ixx=Ixx+(x-cx)*(x-cx)*weight
            Iyy=Iyy+(y-cy)*(y-cy)*weight
            nIxy=nIxy-(x-cx)*(y-cy)*weight       #在惯量矩阵中，反对角线上是-Ixy
        Ixx=Ixx/w_sum; Iyy=Iyy/w_sum; nIxy=nIxy/w_sum
        lambd=0.5*(Ixx+Iyy-np.sqrt((Ixx-Iyy)*(Ixx-Iyy)+4.0*nIxy*nIxy) )
        if np.abs(Iyy) >np.abs(Ixx):
            theta=np.arctan2((lambd-Iyy),nIxy)
        else:
            theta=np.arctan2(nIxy,lambd-Ixx)
        #检查与区域方向的差值，进行纠正
        angle_diff=theta-reg_angle
        while angle_diff<= -np.pi: angle_diff=angle_diff+2*np.pi
        while angle_diff> np.pi: angle_diff=angle_diff-2*np.pi
        if np.abs(angle_diff) >precison: theta=theta+np.pi
        ##计算矩形宽度
        dx=np.cos(theta)
        dy=np.sin(theta)
        l_min=0; l_max=0; w_min=0; w_max=0
        for i in  range(px_count):
            x=region[i,1]; y=region[i,0]
            l= (x-cx)*dx+ (y-cy)*dy     #像素点在矩形主轴上的投影
            w= (x-cx)*(-dy)+(y-cy)*dx   #像素点在幅轴上的投影
            if l>l_max: l_max=l
            if l<l_min: l_min=l
            if w>w_max: w_max=w
            if w<w_min: w_min=w
        rect.x1=cx+l_min*dx     #与主轴正方向相反的端点
        rect.y1=cy+l_min*dy
        rect.x2=cx+l_max*dx     #与主轴正方向相同的端点
        rect.y2=cy+l_max*dy
        rect.cx=cx
        rect.cy=cy
        rect.theta=theta
        rect.dx=dx; rect.dy=dy
        rect.width=w_max-w_min  #注意，这个矩形是非主轴对称的
        rect.precision=precison
        rect.p= p
        return

    def rect_inter_low(self,cur_x,x1,y1,x2,y2):
        """
        给定当前x点，得到此列上与p1->p2线段的交点，
        取交点的y坐标作为下界
        """
        if x1==x2 and y1<y2:     return y1
        if x1==x2 and y1>y2:     return y2
        return y1+(cur_x-x1)*(y2-y1)/(x2-x1)

    def rect_inter_hi(self,cur_x,x1,y1,x2,y2):
        """
        给定当前x点，得到此列上与p1->p2线段的交点，
        取交点的y坐标作为上界
        """
        if x1==x2 and y1<y2: return y2
        if x1==x2 and y1>y2: return y1
        return y1+(cur_x-x1)*(y2-y1)/(x2-x1)

    def get_rect_pxs(self,rect,ll_angles):
        """
        生成并返回给定rect对象包含的像素坐标列表
        """
        px_cnt=0
        alg_px_cnt=0
        height,width=ll_angles.shape[:2]
        ##计算并调整矩形的角点
        dy=np.sin(rect.theta)
        dx=np.cos(rect.theta)
        tcorners=[]
        tcorners.append((rect.y1+dx*rect.width/2, rect.x1-dy*rect.width/2))
        tcorners.append((rect.y2+dx*rect.width/2, rect.x2-dy*rect.width/2))
        tcorners.append((rect.y2-dx*rect.width/2, rect.x2+dy*rect.width/2))
        tcorners.append((rect.y1-dx*rect.width/2, rect.x1+dy*rect.width/2))
        if rect.x1<rect.x2 and rect.y1<=rect.y2 : offset=0
        elif rect.x1>=rect.x2 and rect.y1 <rect.y2: offset=1
        elif rect.x1>rect.x2 and rect.y1>=rect.y2: offset=2
        else: offset =3
        corners=[]
        for i in range(4):
            corners.append(tcorners[(offset+i)%4])
        cur_x=np.ceil(corners[0][1])
        cur_y=np.ceil(corners[0][0])
        ys=-np.inf; ye=-np.inf  #列的起始y坐标，终止y坐标
        while cur_x<corners[2][1] and cur_y>ye:
            ##计算得到当前x列的y坐标下界，上界
            if cur_x < corners[3][1]:
                ys=LSD.rect_inter_low(cur_x, corners[0][1],corners[0][0],corners[3][1],corners[3][0])
            else:
                ys=LSD.rect_inter_low(cur_x, corners[2][1],corners[2][0],corners[3][1],corners[3][0])
            if cur_x <corners[1][1]:
                ye=LSD.rect_inter_hi(cur_x,corners[0][1],corners[0][0],corners[1][1],corners[1][0])
            else:
                ye=LSD.rect_inter_hi(cur_x,corners[1][1],corners[1][0],corners[2][1],corners[2][0])
            ##添加在范围内的像素
            cur_y=np.ceil(ys)
            while cur_y<=ye:
                #是矩形中的点
                if cur_y >= 0 and cur_x >= 0 and cur_y < height and cur_x < width:
                    px_cnt=px_cnt+1
                    if LSD.isaligned(np.int(cur_x),np.int(cur_y), rect.theta,ll_angles,  rect.precision):
                        alg_px_cnt = alg_px_cnt + 1
                cur_y=cur_y+1
            cur_x=cur_x+1
        return px_cnt,alg_px_cnt

    def rect_copy(self,rect_src,rect_dst):
        """
        将src中的信息，拷贝到dst中
        """
        rect_dst.x1=rect_src.x1;
        rect_dst.y1=rect_src.y1    #矩形主轴反方向上的端点
        rect_dst.x2=rect_src.x2
        rect_dst.y2=rect_src.y2    #矩形主轴正方向上的端点
        rect_dst.width=rect_src.width       #宽度
        rect_dst.cx= rect_src.cx     #矩形的中心点
        rect_dst.cy=rect_src.cy
        rect_dst.theta=rect_src.theta        #矩形的方向（弧度）
        rect_dst.precision=rect_src.precision
        rect_dst.p=rect_src.p
        rect_dst.dx=rect_src.dx
        rect_dst.dy=rect_src.dy

    def rect_nfa(self,rect, ll_angles,logNT):
        """
        计算矩形的对数NFA值
        """
        #####计算矩形内总的有效像素个数，以及对齐点的个数
        px_cnt,alg_px_cnt=LSD.get_rect_pxs(rect,ll_angles)
        #####计算NFA
        n=px_cnt; k=alg_px_cnt
        if n==0 or k==0: return -logNT
        if n==k: return -logNT- n*np.log10(rect.p)
        log1Term=LSD.log_gamma(n+1)-LSD.log_gamma(k+1) \
                  -LSD.log_gamma(n-k+1)+k*np.log(rect.p)+(n-k)*np.log(1-rect.p)
        firstTerm = np.exp(log1Term)
        ##在有些情况下计算第一项就足够了
        if np.abs(firstTerm)<2.2e-308:
            if k>n*rect.p:
                return -log1Term/np.log(10)-logNT
            else:
                return -logNT
        ##计算更多项
        bin_tail=firstTerm
        pterm=rect.p/(1-rect.p)
        multiTerm=0.0
        lastTerm=firstTerm
        for i in range(k+1,n+1):
            bin_term=(n-i+1)/i
            multiTerm=bin_term*pterm
            lastTerm=lastTerm*multiTerm
            bin_tail=bin_tail+lastTerm
            if bin_term<1.0:
                err=lastTerm*((1-multiTerm**(n-i+1))/(1-multiTerm)-1)
                if err<0.1*np.abs(-np.log10(bin_tail)-logNT)*bin_tail :
                    break
        return -np.log10(bin_tail)-logNT


    def rect_improve(self,rect,rect_bak,ll_angles, logNT,log_eps):
        """
        当矩形没有意义时，尝试矩形的一些变体来改进NFA值。
        返回值：细化的矩形，以及对数NFA值
        """
        log_nfa=LSD.rect_nfa(rect,ll_angles,logNT);
        if (log_nfa > log_eps):
            return log_nfa

        ##尝试细化精度
        self.rect_copy(rect,rect_bak)
        for n in range(5):
            rect_bak.p=rect_bak.p/2     #精度减小一半
            rect_bak.precision=rect_bak.p*np.pi
            log_nfa_new=self.rect_nfa(rect_bak,ll_angles,logNT)
            if log_nfa_new > log_nfa:
                log_nfa=log_nfa_new
                self.rect_copy(rect_bak,rect)
        if log_nfa > log_eps:
            return log_nfa

        ##尝试细化宽度
        self.rect_copy(rect,rect_bak)
        for n in range(5):
            if (rect_bak.width-0.5)>=0.5:
                rect_bak.width=rect_bak.width-0.5
                log_nfa_new=self.rect_nfa(rect_bak,ll_angles,logNT)
                if log_nfa_new > log_nfa:
                    log_nfa=log_nfa_new
                    self.rect_copy(rect_bak, rect)
        if log_nfa > log_eps:
            return log_nfa

        ##尝试细化一边的宽度
        self.rect_copy(rect,rect_bak)
        for n in range(5):
            if (rect_bak.width-0.5)>=0.5 :
                rect_bak.x1=rect_bak.x1 - rect_bak.dy*0.25
                rect_bak.y1=rect_bak.y1 + rect_bak.dx*0.25
                rect_bak.x2=rect_bak.x2 - rect_bak.dy*0.25
                rect_bak.y2=rect_bak.y2 + rect_bak.dx*0.25
                rect_bak.width = rect_bak.width-0.5
                log_nfa_new=self.rect_nfa(rect_bak,ll_angles,logNT)
                if log_nfa_new>log_nfa:
                    log_nfa=log_nfa_new
                    self.rect_copy(rect_bak, rect)
        if log_nfa >log_eps:
            return log_nfa
        ##尝试细化另一边的宽度
        self.rect_copy(rect,rect_bak)
        for n in range(5):
            if (rect_bak.width-0.5)>=0.5 :
                rect_bak.x1=rect_bak.x1 + rect_bak.dy*0.25
                rect_bak.y1=rect_bak.y1 - rect_bak.dx*0.25
                rect_bak.x2=rect_bak.x2 + rect_bak.dy*0.25
                rect_bak.y2=rect_bak.y2 - rect_bak.dx*0.25
                rect_bak.width = rect_bak.width-0.5
                log_nfa_new=self.rect_nfa(rect_bak,ll_angles,logNT)
                if log_nfa_new>log_nfa:
                    log_nfa=log_nfa_new
                    self.rect_copy(rect_bak, rect)
        if log_nfa >log_eps:
            return log_nfa

        ##尝试再次细化精度
        self.rect_copy(rect,rect_bak)
        for n in range(5):
            rect_bak.p=rect_bak.p/2     #精度减小一半
            rect_bak.precision=rect_bak.p*np.pi
            log_nfa_new=self.rect_nfa(rect_bak,ll_angles,logNT)
            if log_nfa_new > log_nfa:
                log_nfa=log_nfa_new
                self.rect_copy(rect_bak, rect)
        if log_nfa > log_eps:
            return log_nfa

        return log_nfa

    def reduce_region_radius(self,rect,region,px_count,grad_mag,ll_angles,used,reg_angle,precision,p):
        xs=region[0][1]
        ys=region[0][0]
        ##计算区域半径
        rad1=np.sqrt((rect.x1-xs)**2+(rect.y1-ys)**2)
        rad2=np.sqrt((rect.x2-xs)**2+(rect.y2-ys)**2)
        if rad1>rad2: rad=rad1
        else: rad=rad2
        while True:
            rad=rad*0.75
            ##从区域中剔除在半径之外的像素点,并设置为未使用
            i=0
            while i<px_count:
                xx=region[i,1]; yy=region[i,0]
                if np.sqrt((xx-xs)**2+(yy-ys)**2)>rad:
                    used[yy,xx]=0
                    #注意，此时只是用列表尾元素替代当前元素，i值不需要增加
                    region[i]=region[px_count-1]
                    px_count=px_count-1
                else:
                    i=i+1
            if px_count<2: return None,None,False
            self.region2rect(rect,region,px_count,reg_angle,precision,p,grad_mag)
            dist=np.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
            density=px_count/(dist*rect.width)
            if density>=self.density_thresh:
                break
        return px_count,reg_angle,True

    def refine(self,rect,region,px_count,reg_angle,grad_mag,ll_angles,used,precision,p,density_thresh):
        """
        首先判断矩形中对齐点密度,达到阀值标准则判定为有效线段，否则进行角度容差和区域的调整使得符合标准
        返回值：region,reg_angle,rect,passed
        """
        dist=np.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
        density=px_count/(dist*rect.width)
        #如果满足密度标准，不做处理，返回True
        if density >= self.density_thresh:
            return px_count,reg_angle,True
        #####若不满足标准，首先尝试减少角度容差
        ##计算新的平均角度和容差
        xs=region[0,1]     #区域种子（seed）点的坐标
        ys=region[0,0]
        angle_s=ll_angles[ys,xs]
        sum=0.0; s_sum=0.0; n=0
        for i in range(px_count):
            xx=region[i,1]; yy=region[i,0]
            used[yy,xx]=0   #标记为未使用
            dist=np.sqrt((xx-xs)**2+(yy-ys)**2)
            if dist< rect.width:  ##使用seed点在小范围邻域，半径为矩形宽度
                angle=ll_angles[yy,xx]
                ang_d=LSD.angle_diff(angle,angle_s)
                sum=sum+ang_d
                s_sum=s_sum+ang_d**2
                n=n+1
        mean_angle=sum/n
        tau=2.0*np.sqrt((s_sum-2.0*mean_angle*sum)/n+mean_angle**2)     #两倍标准差
        ##重新进行区域生长
        reg_angle,px_count=LSD.region_grow(region[0,1],region[0,0],ll_angles,region,used,tau)
        if px_count<2:
            return None,None,False
        ##重新计算矩形
        LSD.region2rect(rect,region,px_count,reg_angle,precision,p,grad_mag)
        #重新计算密度
        length=np.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
        density=px_count/(length*rect.width)
        if density<density_thresh:     ##如果还达不到标准，则尝试减少seed的邻域半径
            #return self.reduce_region_radius(rect,region,px_count,grad_mag,ll_angles,used,reg_angle,precision,p)
            return LSD.reduce_region_radius(rect, region, px_count, grad_mag,
                                            ll_angles, used, reg_angle, precision, p,density_thresh)
        ##达到这里，说明标准达到了
        return px_count,reg_angle,True

    def LineSegmentDetection(self,image):
        """
        进行线段检测, 返加所有检测到的线段矩形对象列表
        """
        ##减采样
        #scaled_image=self.sub_sample(image)
        scaled_image=LSD.gaussian_sampler(image,self.scale,self.sigma_scale)
        #scaled_image=self.gaussian_sampler(image,self.scale,self.sigma_scale)

        ##计算梯度幅度，角度，梯度伪序
        precision=np.pi*self.ang_thresh/180
        p=self.ang_thresh/180
        rho=self.quant/np.sin(precision)
        #grad_mag,ll_angles,list_p =self.levelline_angle(scaled_image,rho,self.n_bins);
        grad_mag, ll_angles, list_p = LSD.levelline_angle(scaled_image, rho, self.n_bins);
        height,width=ll_angles.shape[:2]
        """
        for i in range(height*width-1,height*width-20,-1):
            print(list_p[i],end='\t')
            print()
        """
        print("list_p=%d"%len(list_p))
        used=np.zeros([height,width],dtype=np.int)
        ##logNT就是NFA公式中与图像宽高有关的项的对数
        logNT = 2.5 * ( np.log10( width ) + np.log10(height) ) +np.log10(11.0)
        min_region_size=np.int(-logNT/np.log10(p) )
        ##查找线段
        reg=np.zeros((height*width,2),dtype=np.int)
        ls_count=0
        rect_list=[]
        for px in list_p:
            #px=tuple(px)
            if used[px[0],px[1]]==1 or ll_angles[px[0],px[1]]==self.NOTDEF:
                continue
            #对当前像素执行区域生长算法
            #reg_angle,px_count=self.region_grow(px[1],px[0],ll_angles,reg,used,precision) #tested
            reg_angle,px_count=LSD.region_grow(px[1],px[0],ll_angles,reg,used, precision)
            if px_count <min_region_size:
                continue
            #为区域构建矩形
            rect=LSD.Rect()
            rect_bak=LSD.Rect()
            #rect=Rect()
            LSD.region2rect(rect,reg,px_count,reg_angle,precision,p,grad_mag)
            #self.region2rect(rect,reg,px_count,reg_angle,precision,p,grad_mag)
            #####refine
            #px_count,reg_angle,passed=self.refine(rect,reg,px_count,reg_angle,grad_mag,ll_angles,used,precision,p,self.density_thresh)
            px_count,reg_angle,passed=LSD.refine(rect,reg,px_count,reg_angle,grad_mag,ll_angles,used,precision,p,self.density_thresh)
            if not passed:
                continue
            #计算NFA值
            log_nfa=LSD.rect_improve(rect,rect_bak,ll_angles,logNT,self.log_eps)
            if log_nfa<=self.log_eps:
                continue
            ls_count=ls_count+1     ##找到并确认新的线段
            ##因为求梯度时是用的2x2掩膜，所以通过偏移校正
            rect.x1=rect.x1+0.5; rect.y1=rect.y1+0.5
            rect.x2=rect.x2+0.5; rect.y2=rect.y2+0.5
            if self.scale!=1.0:
                rect.x1=rect.x1/self.scale; rect.y1=rect.y1/self.scale
                rect.x2=rect.x2/self.scale; rect.y2=rect.y2/self.scale
                rect.width = rect.width/self.scale
            rect_list.append(rect)

        return rect_list